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A SURVEY ON SOME MATHEMATICAL MODELS OF THE 
VERY LOW FREQUENCY WAVE PROPAGATION 

I. INTRODUCTION 


1. Historical Backgroun d 

The development of the theory of the very low frequency wave propagation 
may be divided roughly into two periods. The first period was ushered in with 
the introduction of long distance communications by means of the wireless at 
the turn of the century. The underlying principle in almost all the earlier 
analyses was the diffraction concept advanced by Rayleigh in 1904 [1] . 

In such a diffraction model, the existence of the ionosphere was ignored, 
and the earth was considered as a homogeneous conducting sphere. Among the 
earlier investigators were Poincare [ 2] , MacDonald [ 3] , Nicholson [4] , Som- 
merfeld [ 5] , Love [6] , March [7] , Rybezynski [8] , and others. Most of the 
theoretical investigations made use of the theory of the zonal harmonics, either 
by a direct summation or by forming an integral from the series, to explain the 
facts concerning the propagation of radio waves around the earth's curvature. 
However, the series then formulated was difficult to sum because of its slow 
convergence for the radio problem. 

Poincare [2j and Nicholson [4] attempted to replace the series by an inte- 
gral and then obtained an approximate value for the integral by means of the 
calculus of residues. Their scheme, though substantially sound, was exceedingly 
elaborate, and their approximations became invalid in the vicinity of the antipodes 
of the transmitter. 

The method used by MacDonald [3] was to approximate the original series 
by a modified series which in turn was replaced by an integral. This procedure 
was questioned from a purely mathematical point of view, though worked in most 
practical circumstances. 

Love [ 6 ] undertook the laborious task of evaluating the series for the radio 
problem by numerically adding the terms of the series taken in groups. His 
memoir contained a complete bibliography of the problem up to 1915. 

The results of the aforementioned authors did not explain the propagation 
phenomenon round the earth's curvature; nevertheless, a correct conclusion was 
reached that the diffracted field decayed exponentially from the transmitter. 
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Sommerfeld [ 5 ] investigated the problem using a planar model and con- 
cluded that the field decayed approximately as 1/ j£l or at most 1/d where d is 
the distance of propagation from the transmitter, which checked fairly well with 
the observed facts. Upon his suggestion, March [7] reformulated the problem 
using a spherical model. He cast his solution in an integral in order to obtain 
an answer more consistent with observed facts. Indeed, he found that the field 
decayed approximately as 1/ fr sin 6 ~ l/d where d is the angular distance be- 
tween the transmitter and the receiver. March concluded that the waves were 
propagated without exponential attenuation contrary to the results obtained by 
other investigators. His method was severely questioned. However, von 
Rybezynski [8] extended March’s method and pointed out that the exponential 
factor was neglected by March. 

In reference to these investigations, Smith [9] remarked that "many years 
of work of the most famous mathematicians were inadequate to derive this task, 
which at the first glance seems very simple, from that state with respect to 
which Nicholson said that out of all mathematical problems it is the only one 
that has caused such a great divergence of opinions." It was against this back- 
ground that Watson [10] , in 1918, succeeded in his first of two classical papers 
on VLF wave propagation in making the final conclusion that we could not explain 
the fact of the long distance propagation of long waves around the terrestrial 
sphere using only the diffraction model. 

The second period was initiated in 1919 by the work of Watson [10]. In the 
second of his two papers cn the radio problem, Watson postulated a concentric 
plasma of high conductivity to account for the ionosphere which is mainly re- 
sponsible for the VLF propagation over great distances. Watson succeeded in 
reformulating the integral used by others by means of now the well known Watson 
transformation. Instead of working with the slowly convergent series he was 
able to replace it by a rapidly converging series. He introduced the earth iono- 
sphere waveguide concept. The radio waves excited by a vertical Hertzian dipole 
propagated in the atmosphere between the earth and the ionosphere in a way 
similar to the microwave in modern wave guides. 

Among the earlier investigators in the second period in the development of 
the theory of the propagation of LF and VLF waves, aside from Watson, are 
Rydbeck [ 11] , Van der Pol and Bremmer [12] , Eckersley [13] , Bremmer [14], 
Kendrick [15] , Weyrick [16] , Brekhovskikh [17] , and Al'pert [18] , Fock [19] , 
and others. The results of their analyses indicated needs for further refine- 
ment of the model such as the introduction of terrestrial sphere stratification, 
anisotropy, continuous ionospheric stratification, polarization coupling and 
anisotropy at the ionosphere and the effects of ions, both heavy and light. 
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2. Recent Investigations 


Among the most recent and active investigators are Budden [20], Wait [21] , 
Johler [22], and Kranushkin [23]. 

Wait is a proponent of the mode theory (or the residue theory). To account 
for the behaviors of the fields in the presence of an ionosphere, he considered 
successively the ionosphere as sharply bounded, gradual, and then stratified; 
isotropic and the anisotropic; uniform and non-uniform. The stumbling block of 
the mode theory was the transcendental model equation in Hankel functions of 
complex order for determining the modes of the propagating waves. In its solu- 
tions, various approximations were employed. In view of this difficulty, Johler 
attempted the solution of the VLF propagation problem using only the original 
zonal harmonic series. With the improved summation method, the series could 
readily be summed on a large computer in spite of a large number of terms 
needed. Budden, on the other hand, advocated a direct and numerical solution 
of the full- wave equation taking into account of the coupling caused by the 
anisotropic and nonuniform nature of the ionosphere. 

Kranushkin resolved the problem using the theory of non -self- adjoint dif- 
ferential operator eigenfunction expansion. He called his method the normal 
wave solution. He considered the ionosphere as stratified media with the bound- 
aries being weakly angle dependent. 

We will consider in the next section the four mathematical models currently 
being used; namely, 

i. the full-wave solution, 

ii. the mode theory, 

iii. the zonal harmonic series, 

iv. the normal wave method. 


3. Theoretical Background 

Most of the theoretical considerations of problems involving the propagation 
of radio waves are based on the assumption that the electromagnetic field satis- 
fies Maxwell's equation in the form 

V x F - - i 

VxH= ioie E 
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In the above equations, a time harmonic factor exp(io>t) ij understood; E and H 
are the complex representations of the electric and magnetic field vectors; 
and e are the permeability and permittivity of the material medium. For most 
non-ferromagnetic media, the permeability equals to that of the free space. For 
the ionosphere, the permittivity is a tensor. If e is put equal to 6 q k , where e Q 
is the permittivity of the free space, then k , the dielectric tensor can be written 
as 


0 o\ 


1 0 

I U (U J - Y J ) 

0 1 / 


/ U J - i 2 Y 2 

i n Y U - i m Y J 

-imYU-'tnY 2 ^ 


x -inYU-UY 1 

U 2 - m J Y 2 

iUl-nnY 1 

i 


\ iinYU-UY 1 

-iUU-.nY 1 

U J -n 2 Y 2 ) 


The standard URSI notation is used here. X is the square of the ratio of plasma 
frequency to wave frequency; U = 1 -iZ, where Z is the ratio of the electron col 


lision frequency to the wave frequency; Y (i , m, n) are the components of the 
vector Y referred to some rectangular cartesian coordinates. Y is the mag- 
nitude of the ratio of the gyromagnetic frequency of electrons to the wave fre- 
quency in the direction of the earth's magnetic field. 

The dielectric tensor is Hermitlan. It implies that the ionosphere is an 
anisotropic medium. The changes in the electron density and collision frequency 
and their variations with height contributing to the inhomogeneity to the property 
of the ionosphere of the VLF wave propagation. 

In all the problems considered in the sequel, the primary source is a 
Hertzian vertical dipole. 
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II. THE DIFFRACTION MODEL 




ft 


* 


fc“ 



V 



We will consider first the diffraction problem together with Van der Pol 
and Bremmer. 

An electric dipole is situated at Q, a distance b away from the center of 
the sphere of radius a. The observation point is at P, a distance r away from 
the center of the sphere. The points Q and P are separated by a distance P 

The electromagnetic fields, E and H satisfy MaxwelPs equations. A time 
factor exp (-icut) is used. The analytic problem consists of finding an electric 
Hertzian vector potential function tt r, where r is a position vector, si. oh that 
it satisfies the following differential equations and boundary conditions: 


W (V 2 + k 2 ) tt = 0 (r > a). 
(V 2 + k 2 ) v - 0 ( r < a). 


(2) — (rv) and k 2 v are continuous at r = a, 
o r 


(3) A singularity at Q of the form 


ik, R 


i kj R 


which is further called the 


primary field. 

m 

The electric and magnetic fields are given by 


E = (k 2 + — — V r n) 
r \ 3r 2 / 


. 1 3 * , v 

*77 = — ( r 77 ) 

r 3 r 30 


H, = -IS k 2 — 

* OJ 3 0 


and 


k 2 

K l. 2 


€ 1. 2 2 « 
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The Zonal Harmonic 


The solution of the total Hertzian potential function in the form of the zonal 
harmonic series for r > a is, 


ikj R 


tot ikjR 


GO 

+ y (2n + 1) R tl (l<1 - tj*’ ( k , •») ( k i r> P„ (cos 6)<p. la) 

^ O', a) 


and for r < a is 


»«.»=■£ ^ (2n + l)(R + n S<‘>(k,b)^ B (k 1 r )P n (cos5) 

2 n = 0 


where 




*P„ (x) = ]r {£<'> (X) + Z< J > (x)} 


• W*>-* (-J- fJFH 


(H = Hankel function, J = Bessel function) and where 


“ l°g (*)n .FI - I. 

lx dx n I _ . lx dx 

*- -* at “ k j • ^ 

- ± log {x((» (x))l -fii 

[x dx J, , l_x dx 


log {x\| / 


log {xv (x) } 


j k 3 k „ a 
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It is well known that 


i k, r _ZL, 

=/ (2n + 1) i<>> (k r)* (k b) P (cos 0 ) (r > b) (0.3a) 

i k . R » 

» n = 0 


i k. R 

e 


ik,R 



(2n + 1) (k, b) »!< n (kj r) ? n (cos 0) (r < b) 


(0.3b) 


By substituting 0.3b into 0.1a, 1.1a can be expressed explicitly as a zonal har- 
monic series. Furthermore, on the surface of the earth were r = a, it can be 
reduced to the following 


n, °' ~( k ,°) 3 S <2n 


♦ i) tJ- p b ( cos ^ 


(0.4) 


where 


N (x. y) = 


log {x (x)> 
Lx d x n 


x 3 k j H 


17 ^ log{yw n(y»] 

L J y s kja 

In the reduction to eq. (0.4), use is made of the Wronskian: 


(0.5) 


2 . 


(*) W' n (X) =— ^ 

X 2 

The Watson Transformation and t he Residue Series 

Following Watson, we transform (0.4) into an integral over n leading to 


77 


tot 


1 

(kj a) 3 



n d n 

cos (n 77 ) 


ft 1 -*!/, < k l b > 

< k l a > 


* P , ,, [cos 

M n- 1/2 ' /J 

— 4 


(0.7a) 


where Cj is the contour which enclosed the positive real axis in a counter- 
clockwise direction. By a judicial manipulation of C , eq. (0.7a) is equivalent 
to the following: 
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77 
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i_ [ , " d " _J P [cos (77 - 

Jl COS ("^ {<»>,„ (V, •) N n-1/J 


#)] 


+ 



f (n) d n 


(0.7b) 


where f(n) is the integrand of (0.7a). The line integral in (0.7b) is usually con- 
sidered negligible. L is the contour encloses all the poles in the integrand of 
(0.7a) other than the ones due to cos(n77). 

The first term in (0.7b) enables us to obtain a physical interpretation. Let 
us expand l/cos (n^ ) as 


_1 

cos(nrr ) 


= 2 e 



m = fl 


gifrm (2n + 1) 


We can write (0.7b) over the integration path L as 


( 0 . 8 ) 


where 


CD 



m = 0 


(0.9) 


77 = 

m 


9 «i rr m C (k . b) t 

_ ndn e i7Tn 12 — — 1 

.'l 


(k, a)3 


t.M/i < k l 8 > 


n- 1/2 


[cos (77 - 0)3 


( 0 . 10 ) 


The integral in (0.10) can now be represented as the sum of the residues of 
all the enclosed poles in L, and it is written as 


77 
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4 /7 i e 1 n m 
(k, a) 3 



iTfri (2m *1) 
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( 0 . 11 ) 
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Upon using a well-known asymptotic formula for the spherical harmonics 
of complex order with positive imaginary part, (if 0 is not too near 0 orn ), 


P n- 1/2 ^ COS ^ ^ 


_1 g“in(7r-5)+irr/4 


i2 rr n s i n 9 


( 0 . 12 ) 


7 T 


# 2 yr?, c i7rrn,+3/4 


(k T a) 3 in 0 
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, 1/2 


in s (9 ♦ 2tt m > 
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(0.13) 


The last factor in (0.13) becomes 


in (9 + 2rr m) ( i a -/? ) (9 + 2n m) 
e = e 


(0.14) 


where 


,k s = a s + 1 (fi. > °> 

Hence (0.13) represents a wave travelling round the sphere m times before 
reaching the receiver with a phase velocity and attenuation determined by a s 
and /0 . Because of the fact that a wave travelling m times round the earth 
becomes greatly attenuated, only the m = 0 term contribute significantly to the 
fields. Thus, 


2 yiTr e i3n/4 
(k t a) 3 /sin 0 


ti 0 -!/* < k . b > 
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1/2 


7^0 ^ n 1 - 1/2 


ON 


n -1/2 
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(0.15) 
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3. The Eigenfunction Method 


Van der Pol and Bremmer also consider the possibility of casting the radio 
wave propagation as an eigenvalue problem. In this case eq's (0.1a) and (0.1b) 
are determined by the following considerations: 

(1) (V 2 + k 2 X 2 ) </>j = 0 (r > a) 

(V 2 + k 2 \ 2 -) 0. = 0 (r < a) 


(2) k 2 0. and — ( r 0. ) are continuous at r = a. 
1 B r J 


(3) k 2 f 4>] 

Jr = a 


? d r + k 2 


0? d r = 1 

j 


Jr - 0 


where 0. ’s are the eigenfunctions corresponding to the eigenvalues * s . The 
primary field pertaining to a point source leads to a three dimensional Dirac 
function 4vi3(R)/k 1 at Q. 

The eigenfunctions are of the form 

(A . (\ . k. r) P (cos 0) (r > a) 

B ,'|/ (\ .k r) P (cos 6 1 ) (r<a) 

n , j 1 n ' n , j 1 ' n' /v ' 


0 


n. J 


(0.16) 


where n is a positive integer. The total field satisfying conditions (1) through 
(3) becomes 


od n 


77. 


4 7/1 


tot 
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(0.17a) 


(r > a) 
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where 


A 2 
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(2 n + l) 
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- f-) 10 

x ax/ 


g {x tf 1 ) (x)} 


— (— + - log (y 0 n (y)> 
k 2 \d y 2 y d y / . 


(0.18a) 


and 


B 2 . = 

n. J 
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2 77 k 2 a 3 {0 n (y)} 2 
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K 2 

(— 

1 

±\ 

k 2 
LK ! 

\d x 2 
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(— + i f) 
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(0.18b) 


The difference between the earlier method and the method of the eigenfunc- 
tions is clear. In the latter the summation is done over the integers whereas in 
the former it is done over the complex roots of the modal equation. 


To sum up, we have seen in this section that there are now clearly three 
different methods used in the solution of the VLF propagation problems. They 
are (a) the zonal harmonic series method; (b) the residue series method; and 
(c) the eigenfunction method. In the next two decades, most of the analyses are 
done following these three methods. 


III. THE EARTH -IONOSPHERE MODEL 


To account for the effects on the VLF wave propagation due to the presence 
of the ionosphere, additional boundary conditions were introduced. The iono- 
sphere was first considered as a sharply bounded homogeneous medium, then as 
stratified media. Magneto-ionic theory was also introduced to take care of both 
the isotropic case where the geomagnetic field was ignored and the anisotropic 
case where the geomagnetic field is included. The additional conditions intro- 
duced led to a great deal of complexity; nevertheless, the mathematical models 
employed had their genesis in the diffraction model. 


1. The Mode Theory 

Among the most active proponent of the mode theory is J. R. Wait. W'e shall 
look at the problem with him. 

A. Isotropic, uniform, sharply bounded ionosp here. The ionosphere is 
treated as a concentric homogeneous sharply bounded ionosphere. The source 
is a verticle dipole. The fields are expressed in terms of a Hertzian vector 
which has only a radial component U. In the region between the surface of the 
earth (r = a) and the lower edge of the ionosphere (r = c,ora + h), 


E 

r 



(rll), 


E, 


1 

r 


7)2 


(rU), 



i e u) 



E * = H r - H, = 0. 


( 1 . 1 ) 


The time factor exp (i ^t) is understood. 

e and m a*. e the electric constants and k = {e^i ) 1 2 oj. Subscripts g and i are 
added to the field quantities where reference is made to the region r < a and 
r > c, respectively. To account for the singularity introduced bv the source at 
r = b, the Hertzian potential satisfies the inhomogeneous wave equation, 
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where 


E£ q > = Z ( p H (q > 


at r = c 


(1.9) 
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1 a r q K 


i eo) r J q ( k R r) 


( 1 . 10 ) 


and 


Z( q > = - 


- — [r h( (k. r)] 
1 ar q 


1 ^ ^ r h< 2 ) (k. r) 

q v l 


These can be rewritten upon replacing kr by x as 


( 1 . 11 ) 


— (xU) = i (Z< q ) /t/) U for x = k a 
x a x B 


( 1 . 12 ) 


i a 


— — (xU) = - i (Z( q V77) U for x = k c 
x a x 1 


(1.13) 


The solution for a < r < c is 
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= ■ * k C y (2 q + 1) h< 2 > (k b) h^ 1 ) (k r) 
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(1.17) 


£ 

-5“ ^» h S !) W),n.- i r /7 > 

A Mnxh<*> (X)) x = kc + iZ(« 0 / 7 , 

q ) _ _ ^_2 

i- fexh(') (x)] x , kc + iZ^/v 


(1.18) 


By means of the Watson transformation, the solution becomes 


U = - i k C y \ — h< 2) (kb) (k r) -^-P [cos (tt-8)] (1.19) • 

/ , sin 1 / 7 t v v D v 

1/ 

V 

where All the values in the summand are to be evaluated at the 

poles of f (v) which are the roots of the modal equation 

D v = 0 ( 1 . 20 ) 

If use is made of the Airy integral as expounded by V. A. Fock in place of the 
Hankel function, the radial component of the electric field can be written in the 
form 


- i k a(9 

E r * V 

a ( 9 sin #) 1/2 


apart from a constant factor. And 


V = 



where 


e" ixt [Wj (t - y) + A (t) w 2 (t - y)] 
[wj(t)- q w J(t)][l - A (t) B (t)] 


dt 


x = (k a/2) 1/2 0, y = (2/k a) 1 ' z k (r - a) 


( 1 . 21 ) 


( 1 . 22 ) 


(1.23) 


Now the contour is to enclose the complex poles t n ’s which are solutions of 

1 - A (t) B (t) = 0 (l- 24 ) 
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where 


A (t) = - 


w i ( t - y 0 ) +£ ?i w i - y 0 > 


w' 2 (t - y 0 ) + q. w 2 (t - y 0 ) 


(1.25) 


B(t) = - 


w' 2 (t) -qw 2 (t) 


wi (t) -qw, (t) 


(1.26) 


y 0 = (2/k a) 1/3 k h, = - i (k a/2) 1/3 z/t^, 7j 0 = 120 tt 
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i e n go \‘>4, 
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The rescue or mode series is 
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(1.27) 


It may also be written as 


[wj (t ) - qw (t )] p_A(t)B(t) 

in 1 n 9 t t * t 

_ J r 
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_ 4 (tt x) 1/2 e _i 


IT/ 




(y) (y) a 


n = 0 


(1.28) 


where y = (2/ka) 1/3 kZ h . The functions G are height functions and they are 
normalized to unity for y or y 0 equal to zero. Explicitly 


G n(y> 


f (t ) 

v n, y ' 

f (t n ) 

' n, 0 ' 


(1.29) 


f (t n . y) = Wj (t n - y) + A (t n ) w 2 (t n - y) 


(1.30) 


The function A is a modal excitation factor. It is normalized so that, in the 
limit of zero curvature and perfectly reflecting boundaries, it becomes unity 
for all modes. Explicitly, 


A - V ° 

A "-T 


(t -q*> - 


(tn-yp-q?) K<t n ) -q» 2 (t n )ja ~ 


-1 


(1.31) 
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b. Anisotropic Ionosphere. For the case of an anisotropic ionosphere where 
the earth magnetic field is not ignored, the electromagnetic field components in 
the homogeneous space (aSr^c) are expressed in terms of purely TM and TE 
waves. They are derived from two scalar functions U and V as follows 
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The functions U and V satisfy 


(V 2 + k 2 ) 



= 0 


(1.33) 


in source free region. The source again is a vertical dipole. The primary field 
of such a source is of a purely TM character. It is derived from a single scalar 
function U e satisfying the following inhomogeneous wave equation, 


( VJ + k J ) U = C 8 (r - b) 5 ( £l (1.34) 

2 v r 2 sin 6 
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The boundary conditions, in terms of surface impedance, are 


— (xU) = i A xU 
3 x * 


— (xV) = (i/A ) xV 
Bx * 


x ■ k a 


(1.35) 


where 


and 


>, = z t /7), u =u, + u t , v = W; 


£ = A n i (XV) - i A„ (xU) 


i (xV) = A 21 (xV) - i A 22 (xU) 


(1.36) 


x * k c 


where 


and 


A 11 =Z 11 /7 J’ A l2 =Z n /7 >' A 21 =Z 21 /7 ) 


A 22 = Z 22 /r >' 


The radial field in the concentric waveguide region written in the matrix 
form is 


where 
F 



p -ik «0 

c 

K 

i Cl) I d s 

J>" L 

a (<9 sin 6) l/ 2 

5_ 



(1.37) 


i5u 


-ixt. 


[Wj (t n - y) + A (t n ) w 2 (t n - y)] 


= . 2 (t t x) i/a «" iTr/4 - — 

n«o [wj (t n ) - q w, (t n )] ^ A(t) B(t) 


(1.38) 
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The roots t n 's are determined from the matrix equation 


IA (t)J [B (t)l = I 


(1.39) 


where 


A (t) = 


and 


r ii jJ'i 


R 1 1 R 1 


exp [- i (4 3) (y n - t) 3 2 - i 77 / 2 ] 


W 2 (t) - q w 2 (t) 
wj (t) - q Wj (t) 


0 


B (t) = 


0 


w j (t) - q w 2 (t) 
[w [ (t) - q Wj (t)_ 


(1.40) 


(1.41) 


q = - i (ka '2) 1 3 A , wi th A = 


1 e 0 « 


vl/2 


1 co 


1/2 


O +16 CO 

R R 


1 - 


o -fie aij 

R R 


and 


/cr fie co 

q = - i (k a / 2) l/3 A h , wi th A h = -i i— 

K * ' i e Q co 


1/2 


- 1 


,iRii , R|| » t | Rjl > and ^R, are the reflection coefficients, where 11 and 1 

refer to the direction of the electric vector, relative to the plane of incidence, 
of the incident (first suffix) and reflected (second suffix) wave. 


2. The Zonal Harmonic Series 


The stumbling block of the mode theory is the rather formidable trans- 
cendental modal equation to solve. Recently Wait and Spies made detailed 
computation of the characteristics of the least attenuated modes. However, a 
single mode is not adequate to describe the complete field satisfying Maxwell's 
equations at either intermediate and great distances for daytime models of the 
lower ionosphere. In the case of a nighttime model even more modes are 
needed at all distances. 
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With the speed of the modern digital computers, the argument once held 
against the zonal harmonic solution of the VLF propagation problem is no longer 
convincing. Granted that the zonal harmonic series converges rather slow. 
Johler showed that with improved summation technique, the results of the com- 
putations indicate that full rigor can be achieved with comparative ease at fre- 
quencies less than about 50 kc/s, leaving only the assumed model for the trans- 
mitter and the propagation medium and avoiding the complications of the Watson 
transformation. 


For a simple model of the terrestrial sphere of radius a surrounded by 
concentric ionosphere from r = c to r =°° , the field components are expressed 
in terms of the Hertzian potential functions, n* and v m as follows, with 
exp (iwt) understood, 


r r b sin 9 


3 ( . J 3t7 

— sin 0 — 

B 6 \ B 6 




r b 


B 2 

B r'B 0 


(r v e ) 


H -1- k 2 3 77 * 

^ b fj. Q i co 


E . = 


Uq i 3 77 m 


* b 36 1 


II - 1 3 / . n B 7T n \ 

H = — sin 6 1 

r rb sin 0 B £ \ 3 0/ 


1 B 2 

H„ = — ( r 77 m ) 

r b B r B0 


( 2 . 1 ) 


where 77 c and v m satisfy the wave equations 


77 

(V : f k*7 \ > = 0 

, 77 m 


( 2 . 2 ) 
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The media of propagation are characterized by their electric constants in 
the wave number, k. Thus, for air, with index of refraction, 77 , or dielectric 
constant fc j = 


or 


k 2 = kj 




(2.3) 


For the ground 


k 



O) 


c 


y., . 1 

f O) 


(2.4) 


where cr is the ground conductivity and e is the relative dielectric constant. 
The wave number of the ionosphere is rather intricate. 


k = k a = £’C.: < 2 - 5 > 

There are four distinct propagation components with the complex index of re- 
fraction, rj ™ , where n equals to 0 or e, corresponding to ordinary wave (0), and 
extraordinary wave (e); and m equals to i or r corresponding to upgoing wave 
(i), and downgoing wave (r). They are related to the roots of a Booker type 
quartic in the parameters, which is the result of the simultaneous solution 
of the Langevin equations of charged particles and Maxwell's equation. The 
index of refraction, 17, 

7 ? 2 r £ 2 + sin 2 0. (2.6) 

where <t>. is the angle of incidence on the lower ionosphere plasma. 
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With a vertical electric dipole transmitter, Q77* located in medium of wave num- 
ber kj, the boundary conditions are 

k i < W S + n V ~ k 2 




r ? = w ; 


3 3 

(ttt*) = — (m™) 
dr 1 3 r 2 


(2.7) 


r = a 


k i ( n o + = 1 "Q* 1 “ — ( r77 P 

em r y r 


7 f- r (rn l + rw P =^o 


1 3 

= Q, - — (r 7r!) 

1 me r B r 3 


/Ll n i O) — ( T 7T*) = Q k 2 77® 
0 r 1 Y em 3 3 


( 2 . 8 ) 


where in the quasi- longitudinal approximation, 


0 = -? = ± 


1 k 3 /k . 


F '6> Vk*/kJ - p 


E - sin d 
P = — r = - 

me E* 5 


(2.9) 
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P= sin 2 <£ 2 is the complex direction sine squared of propagation, and the minus 
sign (-) is taken with the ordinary wave, k 3 Q and the plus sign (+) is taken with 
the extraordinary wave k 3 e . The solutions are 


in 

"S=— — r ^ 2n + » (k b) 0 (k rl P n (cos 0), (r < b) 
k 2 r b t—i 


1 n = 0 


ao 

- — (2n + 1) (k j r ) \p n (k l b) P n (cos 0), (r > b) 


k 2 r b 

1 n = 0 


OP 

= = kT 7 b l> n + 1) [b n P „(COS 0 ) 

1 n = 0 


CD 

V 7 S TT-T H (2n + n [b " ? « (k > r) TC " P n < cos V 

k rb 4-1 


1 n = 0 


OP 

77 2 = — — T (2n + 1) a*^' n (k 2 r) P n (cos 0) 
k 2 r b ^ — ■* 


2 n = 0 


w" = 1 > (2n + 11 a m ii fk r) P (cos 6) 

2 UT r u n " 2 


k 2 r b 

2 n = 0 


TT 


00 

= — ~ F (2 n + 1) d* C n (k, r) P (cos ^ 

k 3 rh ‘rrt 


77 = 


UL' 

1 7^ ( 2n + D d : £ n < k 3 r > P „ ( COS ^ 


k 2 r b 
3 n = 0 


(2.10) 


The constants a e , b c and the like are determined by the boundarv conditions. 

n 7 n 

For a stratified ionosphere model, the matching for boundary conditions 
could be a rather formidal task. (For details see Reference 22). 
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3. The Full-wave Theory 


It is seen that when coupling is introduced into the wave propagation problem, 
no matter what mathematical model one uses whether it is based on the residue 
series or the zonal harmonic series, one often resorts to the digital computer 
for numerical computations. Budden is one among many of the English school 
who rely more and more on the direct solution of the full wave equation using 
computers. This approach is formally of complete generality and does not de- 
pend on the assumption of a special ionospheric model. 

In deriving the differential equations, Budden uses the cartesian coordinates 
with the z-axis directed vertically upwards. A plane wave is incident at angle P. 
to the verticle on the ionosphere from below. Let s = sin^, c = cos P., Let the 
ionosphere be stratified such that in each of which the medium is considered 
homogeneous. For the waves in each stratum 3/3 x( ) = - iks ( ); 3/3 y ( ) - 0. 
Then Maxwell’s equations give, 


dE 

y 

d z 


- ik H 

X 



d z 


+ iksE = - ik H 

z y 


-iksE = - i k H 
y * 


(3.1a) 


dH 

y 

d z 



dH 

X 

d z 


ik 


iksH = — D 


- i k s H = D 

y s 2 


(3.1b) 
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From the magnetoionic theory, after the z-components of the fields are 
eliminated, we can write the following matrix equation, 


4— = - i k T e (3.2) 

d z 

where 

f*A 

- F. 

e = y 

H 

x 

\ Hy / (3.3) 

and T is a 4 x 4 matrix whose elements are related to the parameters describ- 
ing the ionosphere, such as the electron density, the electron collision frequency, 
the earth's magnetic field and others. 

The matrix T has four characteristic roots or eigenvalues q A (i = 1, 2, 3, 4) 
which satisfy the characteristic equation 

det (T - ql) = 0 < 3 * 4 ) 

where I is the 4x4 unit matrix. 

In order to solve the matrix equation, it is convenient to introduce a new 
column matrix f, thus 


such that 



e =: S f 


(3.5) 


(3.6) 


where S is some 4x4 matrix whose elements are functions of z to be 
determined. 


25 


Assuming that S is non-singular so that its inverse S" 1 exists. Substitu- 
tion of (3.6) in (3.2) and premultiplication of S _1 gives 

f' + i k S~ 1 TS f = - S' 1 S' f (3.7) 


where the dash indicates the derivative with respect to z. 


The S is so chosen as to make f . the only element of f appearing on the 
left hand side of the i f th equation in (3.7). This means that S -1 TS must be a 
diagonal matrix, this can be done if the characteristic roots q of T are all 
distinct. It can then be shown that 


S" 1 TS 


and eq. (3.7) becomes 


q x 0 0 

o q 2 o 

0 0 q 3 



0 0 


0 



f' + i k A f = - S _1 S' f 


(3.8) 


(3.9) 


Equations (3.9) are the four, first order, coupled differential equations. 
For a homogeneous medium S’ = 0 and (3.9) becomes 


f'+ikAf = 0 (3.10) 

or 

f' + i k q. f. = 0 (i = 1, 2, 3, 4) (3.11) 

The independent characteristic waves in a homogeneous medium are determined 

by 

f; = e lkq ‘ * (i = 1, 2. 3. 4) (3.12) 

Now, premultiply (3.9) by any 4x4 matrix, say F, (3.9) then can be written 
as 

(F f) ' + (i k F A - F') f = - F S' 1 S' f (3.13) 

Choose F so that 

F' = ikFA, (3.14) 
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and 


(F f)' = - F S" 1 S' f 


(3.15) 


and integrate (3.15) in the form 


f = F* 1 - F' 1 



FS ’ 1 S' f dz 


(3.16) 


From (3.14) F may be taken as the diagonal matrix whose diagonal elements 
are exp(ik { q.dz); hence, if f is the fundamental solution of (3.10) namely 


(3.16) becomes 



(3.17) 



J S'f Q dz 


(3.18) 


The evaluation of S _1 S f and its properties can be found in Budden’s book. 
The thing of interest to us here is that we have a set of first order differential 
equations (3.15), which can be integrated numerically choosing various iono- 
spheric models. A work of this kind has been undertaken by Budden and his co- 
workers in the Cavandish Laboratory. Their method of attack is to assume 
some plausible law for the variation of electron density and collision frequency 
with height in the ionosphere, and to calculate the reflecting properties of this 
model ionosphere for long and very long waves. This is repeated for a series 
of models, and those models are selected which have properties most closely 
resembling the experimental results. 


4. The Normal Wave Solution 

In all the previous methods of solution of the VLF wave propagation problem, 
numerical checks were made by assuming a prior knowledge of the parameters 
used. The selections of the parameters though based on sound and scientific 
judgment, were, to say the least, arbitrary. Furthermore, these methods could 
account for only homogeneous tracks (night or day) up to now. The diurnal 
changes in both the phase and the amplitude of the VLF waves remained mostly 
unexplained. In view of the inadequacy in these approaches, Krasnushkin pro- 
posed a mixed method. He suggested to use the data known in the short range 
propagation to evaluate some of the unknown parameters. This is basically an 
inverse method. With the unknown parameters determined, a direct method is 
used to evaluate the fields at great distances. In doing so, some of the diffi- 
culties encountered in evaluating the roots of a transcendental equation as in the 
mode theory could be circumvented. 
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The mathematical models he used are based on the theory of non-self adjoint 
differential operator eigenfunction expansions. He called these functions normal 
waves. 


a. Stratified ionosphere. In a spherically stratified media, let the dielectric 
tensors be as follows, 


r r 


0 

0 


0 

0 

99 

e 9<t> 

— 6^ 

*94> 

€ k 

u 

k = 0, 

1, 2, . 


(4.1) 


where, for the zero layer (k = 0), or (r = 0, r Q ) 


i 4 77-a. 

£ o _ o - e 1 * . i 2 

rr “ 6 99 ~ € + 1 

CO 


f 0 — 0 

’ ~ u 


(4.2) 


for the first layer ( k = 1); (r Q = a, ^ ) 


f l =1 f ^ — o 

rr “ *09 - 1 e 94> ~ 


(4.3) 


and for the other n - 1 layers (k = 2, 3, . . . , N) 


(r,, r 2 ); fiy r 3 ); .... ; (r N _,. r N = 
the tensor l| | is as follows, 


1 CO 


2 


£ rr= 1 + 


a) (v eff - i co) 


99 


1 + 


1 rj} l - 1 ^ 

W [' V ef f “ ‘ "7* + "h-I 


1 “« 6: h 




^ ^ V eff - i ^) 2 + 


(4.4) 
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where 



4 v N e 1 2 
m 


e H, 


O'h — 


m c 


(4.5) 


The components of the field in each layer are 


t-h (sin 6 Hi) = - — e* E* t 22 ii. 


4 77 


r sin 6 B# 


rr r 


1 

r 


— (rift) 

Br * 


— [€ k 
c ve et 


E 





1 

r sin 6 


h ( sin9 ** ) = 



k 

r’ 



_B 
B r 




1 

r 





(4.6) 


They are subject to the following conditions: 

(1) Tangential components are continuous on the spher- 

ical surfaces r k = const., k = 1, 2, . . . , N. 
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(2) In the center of the earth (r = r M =0) the fields are bounded and 

(3) At infinity they tend to zero. 

Now, introduce the Hertzian potential function 


B (r, 0 ) 
A (r, 0 ) 


(4.7) 


where 



1 BB = 1 9A 

r a £?’ ^ r 3 e 


E 

r 


i 

k Q e rf r 2 sin 0 


3_ 

3£ 


si n 


3B\ 

de) 


i 32B 6 e<t> 3a 

k 0 € ee r 9r 30 " re ee de ' 


H = - i 3 / . „ 3 A\ 

r — sin 9 — . 

k Q r 2 sin 0 d& \ 3 QJ 


u _ i 3 2 A 

rl/j - _ • 

k Q r dr 35 

where k 0 = u/c is the wave number of the free space. 
Eq. (4.6) can be expressed in terms of A and B 


i<y) 




(4.8) 


(4.9) 
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Now decompose (4.7) into the product 


B 


Y(r) 

A 


Z ( r) 


0 ( 0 ) 


(4.14) 


Y 

Then i[ k) will act only on | ! and 4™ only on 9 . After replacing 3/3 r by 

d/dr one obtains the differential operator L r defined by ^ k > . Similarly, one 
obtains L Q defined by ^ k ) . Then one has a single operational equation: 



(4.15) 


To find the solution of (4.15), we will find first the solution of the homoge- 
nous equation 


L 

B 

+ K 

3 

r 

A 

9 

A 


= 0 


(4.16) 


From (4.14), we have 


L 


r 


Y 

Z 


Y(r) 

Z(r) 

(r) 

(r) 


L e 0 (9) 

0 ( 0 ) 


or 


L 


Y(r) 

Z(r) 


= x 


Y(r) 

Z(r) 


(4.17) * 


L 9 0 (9) + x 0 (9) = 0 


(4.18) 


The solution of (4.17) exists for eigenvalues x *s (j = 1, 2, . . .) correspond- 
ing to the eigenfunctions 
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Equation (4.18) is f he usual Legendre equation 


where 


1 

si n 0 


_d_ 

A0 


s in 6 



+ l/ i 



+ 1 ) ^ = 0 


(4.19) 


v \ + ^ = x j (4.20) 

For the case of a verticle dipole source placed at r = b we may write the 
solution of (4.15) as follows, 


B 
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np_ y | Y i (O 

cfc>2 jri l z i ^ 


Y j (b) 


N. sin i/. v 
j j 



[cos (77 - O')] 


where P is the dipole moment of the source and 


(4.21) 


N 

) 



d r + f 
‘'o 



(4.22) 


Using the asymptotical approximation of the Legendre function, and expanding 
1/sin v 77 , we write (4.21) as follows, 


00 



Y. (r) If, (b) v\'* 

Z, (r) N 

J J 



exp j\ (n 


+ 1) 277 v . 



where v • = v . + 1/2. The separate terms in the braces of (4.23) denote the sig- 
nals reaching the receiver after passing round the earth n times. Because of the 
attenuation only the waves with n = 0 arriving at the receiver along the shortest 
path should be considered. Thus, 




+ exp i 2 77 n v . + i 
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(4.24) 
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The field components are 
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For calculating the field in the atmospheric layer, a further simplification 
of formula (4.2) is needed. For the details involved, see reference [ 23] . 

b. Stratified medium with relief. To account for the diurnal variations of 
the field, Krasnushkin proposed another model to take into consideration of the 
slow variations of the dielectric tensor (4.1) along the path of radio waves. He 
assumed that 


|| e || = f (r, /x6\ (4.26) 

where 90° - 9 and 0 are latitude and longitude and /z is some small parameter. 
The operational equation now takes the following form 
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+ L 9<fi 
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This equation is not separable due to dependence of L r on 9 and 0. However, 
for small ix one may expand L in powers of ^ . The first approximation yields 
the following eigenvalue equation: 
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From whence one obtains the set of eigenvalues x. (ix9 t /z0),j=l, 2, ... 00 
corresponding to the eigenfunctions 
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For (6 t (p) one obtains the equation 
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If one seeks the solution in the form 

4>j = Uj (9, 0) exp [i S. (@> <P)] (4.32) 

where U j is a slowly varying function of 8 and 0. After substituting (4.32) into 
(4.31), one finally obtains for | | the following expressions: 
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m 

where normalization factors N. (O) and N. (P) are evaluated at the sender and the 
receiver points. The path of the wave is defined by the eikonal function evalu- 
ated from the equation 


and 



If the geodetic line is followed, then 



(4.34) 
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After substituting (4.36) into (4.33), one finally obtains the radial component of 
the electric field 
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IV. SUMMARY 


So far we have seen four mathematical models in the study of the VLF wave 
propagation; namely (1) waveguide mode theory, (2) zonal harmonic series, 

(3) eigenfunction method, and (4) full wave theory. Of the four, the last stresses 
more on the computational aspects of the problem whereas the first three give 
analytic expressions as solutions of the problem. Models one and three are 
more elegant in form, each is capable of explaining the propagation problem in 
simple and tractable terms using mostly one or two modes. Model (1) is based 
on the concept of the classical residue series, whereas model (2) is formulated 
on the modern theory of the non-self adjoint operators. The stumbling block in 
the mode theory is the rather formidable transcendental equation whose solution 
is necessary to account for the various modes needed in evaluating the electro- 
magnetic fields. There is a corresponding transcendental equation in the eigen- 
value method, however, Krasnushkin showed that the difficulty can be circum- 
vented by means cf what he called the mixed method; whereby he used the data 
of the short range fields and some of the known parameters to determine the 
unknown parameters. Then he used the known and the determined parameters 
to evaluate the long range propagation phenomenon. 

In contrast to the mode theory model and the eigenfunction model, the zonal 
series model relies only on the original harmonic series for its solution. The 
speed of the electronic computers offsets the objection of the direct summation 
of the slowly convergent harmonic series. A complete solution is thus made 
possible leaving only the assumed model for the transmitter and the propagation 
medium. The advantage of a descriptive discussion that a simple formula affords 
is lost. Perhaps when all the details of the attributes of the ionosphere are in- 
cluded, graphical solutions based on numerical computation are the only answers. 
In which case, the full wave theory model as is handled by Budden may be the 
best approach. Unfortunately, Budden's is a planar model. To satisfy the long 
range propagation, the earth's curvature must be taken into consideration. 

The abundant experimental materials cumulated over the past fifty years in 
the field of VLF propagation have established many facts that the current theories 
are still incapable of explaining, such as 

(1) the diurnal and the seasonal variations of the long distance VLF field; 

(2) the dependence of both the field intensity and the phase on a hetero- 
genous path of various illuminations; 

(3) the effect of the solar flare on both the field strength and the phase. 

It appears that some work could be done on the formulation of a new model 
which does not depend only harmonically on time in order to resolve some of the 
aforementioned observed facts. 
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